{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from models.BootstrapDY import estimate_var, make_dynamics, make_var_mbb, make_boot_dynamics, rolling_window, static, DieboldYilmaz2012\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "data = pd.read_csv('dy2012.csv', parse_dates=['Index'], index_col=0)\n",
    "#data = data.loc['2007-1-1':'2008-1-1']\n",
    "YY = data.values"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "nlags = 4\n",
    "pp = nlags\n",
    "Horizon = 10\n",
    "TT = YY.shape[0] - nlags\n",
    "forecast_horizon = Horizon"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "outputs": [],
   "source": [
    "from statsmodels.tsa.api import VAR\n",
    "model = VAR(YY)\n",
    "results = model.fit(nlags)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "statsA = results.coefs\n",
    "for i in range(4):\n",
    "\tstatsA[i] = statsA[i].T\n",
    "statsA = statsA.reshape(16,4)\n",
    "statsA = np.vstack((results.intercept.reshape((1,4)), statsA))\n",
    "\n",
    "statsU = results.resid\n",
    "stats_sigmacomp = results.sigma_u\n",
    "statscovUU = (statsU.T @ statsU) / TT\n",
    "\n",
    "statsPhi = results.ma_rep(9)\n",
    "for I in range(10):\n",
    "\tstatsPhi[i] = statsPhi[i].T\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [],
   "source": [
    "numObs, KK = YY.shape\n",
    "fTT = numObs - nlags\n",
    "\n",
    "fyy = YY[nlags:, :]\n",
    "fxx = np.zeros((fTT, 1 + nlags*KK))\n",
    "fxx[:, 0] = 1\n",
    "for t in range(fTT):\n",
    "    for p in range(nlags):\n",
    "        fxx[t, 1+p*KK:1+KK+p*KK] = YY[t+nlags-p-1, :]\n",
    "A_est, U_est, covUU_est = estimate_var(fyy, fxx)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 进行系数的比较(使用sovle版)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[9.21263066e-13, 1.32560629e-13, 7.62057084e-13, 7.90478794e-13],\n       [4.55191440e-15, 9.23898864e-02, 1.91220139e-02, 4.26735851e-02],\n       [9.23898864e-02, 5.35682609e-15, 4.56085347e-02, 2.36257513e-02],\n       [1.91220139e-02, 4.56085347e-02, 4.74620343e-15, 2.31236051e-02],\n       [4.26735851e-02, 2.36257513e-02, 2.31236051e-02, 2.88935542e-14],\n       [4.88498131e-14, 2.27184873e-02, 2.48718457e-03, 3.18704010e-02],\n       [2.27184873e-02, 3.09474668e-14, 3.80694633e-02, 3.46951451e-02],\n       [2.48718457e-03, 3.80694633e-02, 1.28785871e-14, 3.35262082e-02],\n       [3.18704010e-02, 3.46951451e-02, 3.35262082e-02, 4.02455846e-14],\n       [4.62685446e-14, 5.27943399e-03, 3.40257466e-02, 1.32861261e-02],\n       [5.27943399e-03, 1.76525461e-14, 4.77290233e-02, 2.59570507e-02],\n       [3.40257466e-02, 4.77290233e-02, 2.60902411e-15, 3.35859918e-02],\n       [1.32861261e-02, 2.59570507e-02, 3.35859918e-02, 3.13360449e-14],\n       [2.21767049e-14, 5.94488533e-02, 2.72907220e-02, 1.47860703e-02],\n       [5.94488533e-02, 7.85482790e-15, 2.11226387e-02, 7.82172730e-02],\n       [2.72907220e-02, 2.11226387e-02, 2.98094882e-14, 2.20354062e-02],\n       [1.47860703e-02, 7.82172730e-02, 2.20354062e-02, 1.81382687e-14]])"
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(A_est - statsA) #A有时差距较小有时差距较大"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[1.33226763e-13, 1.65201186e-13, 2.48689958e-14, 1.65201186e-13],\n       [3.90798505e-14, 1.06581410e-13, 9.05941988e-14, 3.01980663e-14],\n       [3.55271368e-15, 1.61648472e-13, 9.59232693e-14, 3.37507799e-14],\n       ...,\n       [1.06581410e-14, 2.66453526e-14, 1.17239551e-13, 2.84217094e-14],\n       [2.84217094e-14, 7.28306304e-14, 2.66453526e-14, 5.32907052e-15],\n       [1.42108547e-14, 6.39488462e-14, 3.55271368e-14, 8.17124146e-14]])"
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(U_est - statsU) #扰动项差距较小"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[1.11022302e-16, 5.55111512e-17, 1.38777878e-17, 8.32667268e-17],\n       [5.55111512e-17, 2.22044605e-16, 5.55111512e-17, 5.55111512e-17],\n       [1.38777878e-17, 5.55111512e-17, 2.22044605e-16, 2.77555756e-17],\n       [8.32667268e-17, 5.55111512e-17, 2.77555756e-17, 0.00000000e+00]])"
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(covUU_est - statscovUU) # 如果按公式计算协方差较小"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0.0042642 , 0.00122055, 0.00046343, 0.00076953],\n       [0.00122055, 0.00545408, 0.00041039, 0.00116751],\n       [0.00046343, 0.00041039, 0.00876819, 0.00039462],\n       [0.00076953, 0.00116751, 0.00039462, 0.00443586]])"
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(covUU_est - stats_sigmacomp) #如果直接拿到他给的协方差那差距较大"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0.0042642 , 0.00122055, 0.00046343, 0.00076953],\n       [0.00122055, 0.00545408, 0.00041039, 0.00116751],\n       [0.00046343, 0.00041039, 0.00876819, 0.00039462],\n       [0.00076953, 0.00116751, 0.00039462, 0.00443586]])"
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(statscovUU - stats_sigmacomp) #dy的程序是使用的sigmacomp"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "$\\frac{1}{T - Kp - 1} Y^\\prime (I_T - Z (Z^\\prime Z)^{-1} Z^\\prime) Y$"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.]])"
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.sigma_u - ((statsU.T @ statsU) / (TT - KK*nlags -1)) #statsmodel这里进行了自由度调整"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[-1.11022302e-16, -5.55111512e-17,  1.38777878e-17,\n        -5.55111512e-17],\n       [-5.55111512e-17, -2.22044605e-16, -4.16333634e-17,\n         5.55111512e-17],\n       [ 1.38777878e-17, -4.16333634e-17,  0.00000000e+00,\n        -4.16333634e-17],\n       [-5.55111512e-17,  5.55111512e-17, -4.16333634e-17,\n         1.11022302e-16]])"
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.sigma_u - ((U_est.T @ U_est) / (TT - KK*nlags -1))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 进行系数的比较使用statsmodel的版本"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.]])"
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(A_est - statsA) #A有时差距较小有时差距较大"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       ...,\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.]])"
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(U_est - statsU) #扰动项差距较小"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0.00000000e+00, 0.00000000e+00, 1.38777878e-17, 5.55111512e-17],\n       [0.00000000e+00, 2.22044605e-16, 2.77555756e-17, 5.55111512e-17],\n       [1.38777878e-17, 2.77555756e-17, 2.22044605e-16, 1.38777878e-17],\n       [5.55111512e-17, 5.55111512e-17, 1.38777878e-17, 1.11022302e-16]])"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(covUU_est - stats_sigmacomp) #如果直接拿到他给的协方差那差距较大"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.],\n       [0., 0., 0., 0.]])"
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.sigma_u - ((U_est.T @ U_est) / (TT - 17))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 直接利用A和U进行测试"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0.8772493 , 0.08150053, 0.00489427, 0.03635589],\n       [0.09122472, 0.82270893, 0.03492635, 0.05114001],\n       [0.00331561, 0.02960581, 0.95160661, 0.01547196],\n       [0.05509739, 0.0769218 , 0.02148512, 0.84649569]])"
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "make_dynamics(A_est, covUU_est, 4, 10) # 计算结果仍然有差距"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 继续测试Phi的计算"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "outputs": [],
   "source": [
    "AA = A_est\n",
    "\n",
    "NN, KK = AA.shape                                       # AA这里是行观点系数\n",
    "#coefs = AA[1:, :].reshape((pp, KK, KK))                # 1.提取非截距项部分 2.为了计算方便,将系数矩阵设为3维\n",
    "\n",
    "coefs = np.zeros((pp, KK, KK))\n",
    "for i in range(pp):\n",
    "    coefs[i, :, :] = AA[1+i*KK:1+(i+1)*KK, :]\n",
    "\n",
    "Phis = np.zeros((Horizon, KK, KK))\n",
    "Phis[0] = np.eye(KK)\n",
    "\n",
    "# recursively compute Phi matrices\n",
    "for i in range(1, Horizon):\n",
    "    for j in range(1, i+1):\n",
    "        if j > pp:\n",
    "            break\n",
    "\n",
    "        Phis[i] += coefs[j-1] @ Phis[i-j]  # 注意这里Phis的下标没有位移(从Phi0开始), 但是AA的下标有位移"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "outputs": [],
   "source": [
    "statsPhi = results.ma_rep(9)\n",
    "for i in range(10):\n",
    "\tstatsPhi[i] = statsPhi[i].T\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "outputs": [
    {
     "data": {
      "text/plain": "0.0"
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(Phis - statsPhi).sum()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Phi也没有问题"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "outputs": [],
   "source": [
    "AA = A_est\n",
    "\n",
    "NN, KK = AA.shape                                       # AA这里是行观点系数\n",
    "#coefs = AA[1:, :].reshape((pp, KK, KK))                # 1.提取非截距项部分 2.为了计算方便,将系数矩阵设为3维\n",
    "\n",
    "coefs = np.zeros((pp, KK, KK))\n",
    "for i in range(pp):\n",
    "    coefs[i, :, :] = AA[1+i*KK:1+(i+1)*KK, :]\n",
    "\n",
    "Phis = np.zeros((Horizon, KK, KK))\n",
    "Phis[0] = np.eye(KK)\n",
    "\n",
    "# recursively compute Phi matrices\n",
    "for i in range(1, Horizon):\n",
    "    for j in range(1, i+1):\n",
    "        if j > pp:\n",
    "            break\n",
    "\n",
    "        Phis[i] += coefs[j-1] @ Phis[i-j]  # 注意这里Phis的下标没有位移(从Phi0开始), 但是AA的下标有位移\n",
    "\n",
    "#================================================================================\n",
    "CovUU = covUU_est\n",
    "deviations = np.sqrt(np.diag(CovUU))# 这个没有问题\n",
    "\n",
    "# for i in range(10):\n",
    "# \tPhis[i] = Phis[i].T\n",
    "#--------------\n",
    "numerators = np.zeros((Horizon, KK, KK))\n",
    "denominator = np.zeros((Horizon, KK, KK))\n",
    "for h in range(Horizon):\n",
    "    numerators[h] = Phis[h].T @ CovUU\n",
    "    denominator[h] = numerators[h] @ Phis[h]\n",
    "    numerators[h] = np.square(numerators[h]/ deviations)\n",
    "\n",
    "DD = np.zeros((KK, KK))\n",
    "for i in range(KK):\n",
    "    for j in range(KK):\n",
    "        DD[i, j] = np.sum(numerators[:, i, j])/np.sum(denominator[:, i, i])\n",
    "#---------------\n",
    "\n",
    "\n",
    "DD_tilde = DD/DD.sum(axis=1).reshape(KK,-1)*100"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "outputs": [],
   "source": [
    "sigma_u = np.asarray(results.sigma_u)\n",
    "sd_u = np.sqrt(np.diag(sigma_u))\n",
    "\n",
    "\n",
    "#-----------------\n",
    "fevd = results.fevd(forecast_horizon, sigma_u / sd_u)\n",
    "fe = fevd.decomp[:, -1, :]\n",
    "#-----------------\n",
    "\n",
    "\n",
    "fevd = (fe / fe.sum(1)[:, None] * 100)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[88.75700164,  7.29118463,  0.34532792,  3.6064858 ],\n       [10.21354465, 81.44571151,  2.72697368,  5.61377016],\n       [ 0.46811796,  3.69595287, 93.69418926,  2.14173991],\n       [ 5.69157879,  7.02601676,  1.54775915, 85.7346453 ]])"
     },
     "execution_count": 176,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fevd"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[1.42108547e-14, 8.88178420e-16, 3.33066907e-16, 3.10862447e-15],\n       [0.00000000e+00, 1.42108547e-14, 8.88178420e-16, 2.66453526e-15],\n       [5.55111512e-17, 0.00000000e+00, 0.00000000e+00, 4.44089210e-16],\n       [8.88178420e-16, 6.21724894e-15, 0.00000000e+00, 1.42108547e-14]])"
     },
     "execution_count": 177,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fabs(DD_tilde - fevd)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "outputs": [],
   "source": [
    "def lrange(*args, **kwargs):\n",
    "    return list(range(*args, **kwargs))\n",
    "\n",
    "irfobj = results.irf(var_decomp=(sigma_u / sd_u), periods=10)\n",
    "orth_irfs = irfobj.orth_irfs\n",
    "\n",
    "# cumulative impulse responses\n",
    "irfs = (orth_irfs[:10] ** 2).cumsum(axis=0)\n",
    "\n",
    "neqs = results.neqs\n",
    "rng = lrange(neqs)\n",
    "mse = results.mse(10)[:, rng, rng]\n",
    "\n",
    "# lag x equation x component\n",
    "fevd = np.empty_like(irfs)\n",
    "\n",
    "for i in range(10):\n",
    "    fevd[i] = (irfs[i].T / mse[i]).T\n",
    "\n",
    "# switch to equation x lag x component\n",
    "#decomp = fevd.swapaxes(0, 1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[0.69406084, 0.19866282, 0.07542932, 0.12525187],\n        [0.19866282, 0.88773112, 0.06679637, 0.19003002],\n        [0.07542932, 0.06679637, 1.42715106, 0.06422977],\n        [0.12525187, 0.19003002, 0.06422977, 0.72200096]],\n\n       [[0.71726523, 0.21086386, 0.07762947, 0.13556694],\n        [0.21086386, 0.92522197, 0.09989765, 0.19786279],\n        [0.07762947, 0.09989765, 1.48927624, 0.07490291],\n        [0.13556694, 0.19786279, 0.07490291, 0.73020439]],\n\n       [[0.78261253, 0.23370429, 0.07995612, 0.15431813],\n        [0.23370429, 0.97427757, 0.12628964, 0.2116062 ],\n        [0.07995612, 0.12628964, 1.5816593 , 0.09332389],\n        [0.15431813, 0.2116062 , 0.09332389, 0.76725946]],\n\n       [[0.85317743, 0.2698687 , 0.06926143, 0.17331434],\n        [0.2698687 , 1.02196557, 0.13101231, 0.23299451],\n        [0.06926143, 0.13101231, 1.69132926, 0.10093577],\n        [0.17331434, 0.23299451, 0.10093577, 0.80573988]],\n\n       [[0.94045189, 0.30923633, 0.07197224, 0.1972989 ],\n        [0.30923633, 1.11821147, 0.19374696, 0.26337299],\n        [0.07197224, 0.19374696, 1.83387395, 0.13730869],\n        [0.1972989 , 0.26337299, 0.13730869, 0.84210924]],\n\n       [[0.98444108, 0.33271292, 0.06952747, 0.2118588 ],\n        [0.33271292, 1.15520034, 0.22149974, 0.27822534],\n        [0.06952747, 0.22149974, 1.90125637, 0.15162861],\n        [0.2118588 , 0.27822534, 0.15162861, 0.85656501]],\n\n       [[1.03418221, 0.35740251, 0.06540288, 0.22665428],\n        [0.35740251, 1.19070286, 0.24605777, 0.2937144 ],\n        [0.06540288, 0.24605777, 1.96600415, 0.16601709],\n        [0.22665428, 0.2937144 , 0.16601709, 0.87223308]],\n\n       [[1.07844559, 0.38103029, 0.06125979, 0.23996286],\n        [0.38103029, 1.22009087, 0.26634066, 0.30774115],\n        [0.06125979, 0.26634066, 2.02290186, 0.17944297],\n        [0.23996286, 0.30774115, 0.17944297, 0.88457293]],\n\n       [[1.11781577, 0.40192754, 0.059341  , 0.25223498],\n        [0.40192754, 1.25087273, 0.29181306, 0.32196036],\n        [0.059341  , 0.29181306, 2.0759561 , 0.19364168],\n        [0.25223498, 0.32196036, 0.19364168, 0.89508382]],\n\n       [[1.15067023, 0.41998803, 0.0566537 , 0.26267098],\n        [0.41998803, 1.273563  , 0.30987657, 0.3332263 ],\n        [0.0566537 , 0.30987657, 2.1170029 , 0.20392287],\n        [0.26267098, 0.3332263 , 0.20392287, 0.90308569]]])"
     },
     "execution_count": 152,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "steps = 10\n",
    "ma_coefs = results.ma_rep(9)\n",
    "\n",
    "k = len(results.sigma_u)\n",
    "forc_covs = np.zeros((steps, k, k))\n",
    "\n",
    "prior = np.zeros((k, k))\n",
    "for h in range(steps):\n",
    "    # Sigma(h) = Sigma(h-1) + Phi Sig_u Phi'\n",
    "    phi = ma_coefs[h]\n",
    "    var = phi @ results.sigma_u @ phi.T\n",
    "    forc_covs[h] = prior = prior + var\n",
    "forc_covs"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[ 1.        ,  0.        ,  0.        ,  0.        ],\n        [ 0.        ,  1.        ,  0.        ,  0.        ],\n        [ 0.        ,  0.        ,  1.        ,  0.        ],\n        [ 0.        ,  0.        ,  0.        ,  1.        ]],\n\n       [[ 0.18350846, -0.02424872, -0.00363372,  0.02199887],\n        [ 0.06814117,  0.16377225,  0.05259282,  0.00431903],\n        [-0.02275573,  0.09820135,  0.18333308,  0.04829864],\n        [ 0.06467245, -0.01930672,  0.02517504,  0.06996402]],\n\n       [[ 0.30767075, -0.00641004, -0.01285889,  0.00932753],\n        [ 0.04825777,  0.2049666 ,  0.04305291,  0.01901118],\n        [-0.01182419,  0.0205726 ,  0.24748445,  0.05634644],\n        [ 0.05040629, -0.02434135,  0.01578916,  0.21592331]],\n\n       [[ 0.30996154,  0.02198823, -0.022486  ,  0.01239301],\n        [ 0.08674005,  0.19417692, -0.00985863,  0.02298701],\n        [-0.05701899,  0.04514848,  0.27434849, -0.00235021],\n        [ 0.02811824,  0.02500114,  0.00771659,  0.21441709]],\n\n       [[ 0.34249729,  0.03101059, -0.01540619,  0.01247897],\n        [ 0.04125081,  0.30600514,  0.06310066, -0.01388488],\n        [-0.05830259,  0.10976013,  0.29132875,  0.07670475],\n        [ 0.03594038,  0.04454557,  0.03331973,  0.18771237]],\n\n       [[ 0.24510575,  0.01204579, -0.0181521 ,  0.01619914],\n        [ 0.07675677,  0.16014394,  0.04204776,  0.01400784],\n        [-0.05022241,  0.09040741,  0.19818061,  0.0417974 ],\n        [ 0.05180726,  0.01586668,  0.02372821,  0.10911875]],\n\n       [[ 0.25751386,  0.02322018, -0.02223925,  0.01355798],\n        [ 0.06940312,  0.16041235,  0.03957897,  0.01481425],\n        [-0.05329887,  0.07512706,  0.19710291,  0.04718154],\n        [ 0.04302939,  0.02578756,  0.02140216,  0.11633482]],\n\n       [[ 0.23994421,  0.02901359, -0.02240226,  0.01257221],\n        [ 0.07576793,  0.14245516,  0.03192275,  0.00810255],\n        [-0.06011348,  0.08499565,  0.18036203,  0.04392574],\n        [ 0.04035753,  0.035455  ,  0.02258904,  0.09202785]],\n\n       [[ 0.2251054 ,  0.03012719, -0.01962293,  0.01209812],\n        [ 0.06344208,  0.14752444,  0.04596839,  0.00708544],\n        [-0.05926704,  0.09756479,  0.16733417,  0.04848414],\n        [ 0.04094871,  0.03783219,  0.02646091,  0.0750044 ]],\n\n       [[ 0.20561136,  0.02634501, -0.02096767,  0.01240291],\n        [ 0.07068361,  0.11640777,  0.03766458,  0.01207825],\n        [-0.05526942,  0.08698693,  0.14711928,  0.0407007 ],\n        [ 0.04214355,  0.03201595,  0.02239741,  0.06205465]]])"
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ma_coefs"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[ 1.        ,  0.        ,  0.        ,  0.        ],\n        [ 0.        ,  1.        ,  0.        ,  0.        ],\n        [ 0.        ,  0.        ,  1.        ,  0.        ],\n        [ 0.        ,  0.        ,  0.        ,  1.        ]],\n\n       [[ 0.18350846,  0.06814117, -0.02275573,  0.06467245],\n        [-0.02424872,  0.16377225,  0.09820135, -0.01930672],\n        [-0.00363372,  0.05259282,  0.18333308,  0.02517504],\n        [ 0.02199887,  0.00431903,  0.04829864,  0.06996402]],\n\n       [[ 0.30767075,  0.04825777, -0.01182419,  0.05040629],\n        [-0.00641004,  0.2049666 ,  0.0205726 , -0.02434135],\n        [-0.01285889,  0.04305291,  0.24748445,  0.01578916],\n        [ 0.00932753,  0.01901118,  0.05634644,  0.21592331]],\n\n       [[ 0.30996154,  0.08674005, -0.05701899,  0.02811824],\n        [ 0.02198823,  0.19417692,  0.04514848,  0.02500114],\n        [-0.022486  , -0.00985863,  0.27434849,  0.00771659],\n        [ 0.01239301,  0.02298701, -0.00235021,  0.21441709]],\n\n       [[ 0.34249729,  0.04125081, -0.05830259,  0.03594038],\n        [ 0.03101059,  0.30600514,  0.10976013,  0.04454557],\n        [-0.01540619,  0.06310066,  0.29132875,  0.03331973],\n        [ 0.01247897, -0.01388488,  0.07670475,  0.18771237]],\n\n       [[ 0.24510575,  0.07675677, -0.05022241,  0.05180726],\n        [ 0.01204579,  0.16014394,  0.09040741,  0.01586668],\n        [-0.0181521 ,  0.04204776,  0.19818061,  0.02372821],\n        [ 0.01619914,  0.01400784,  0.0417974 ,  0.10911875]],\n\n       [[ 0.25751386,  0.06940312, -0.05329887,  0.04302939],\n        [ 0.02322018,  0.16041235,  0.07512706,  0.02578756],\n        [-0.02223925,  0.03957897,  0.19710291,  0.02140216],\n        [ 0.01355798,  0.01481425,  0.04718154,  0.11633482]],\n\n       [[ 0.23994421,  0.07576793, -0.06011348,  0.04035753],\n        [ 0.02901359,  0.14245516,  0.08499565,  0.035455  ],\n        [-0.02240226,  0.03192275,  0.18036203,  0.02258904],\n        [ 0.01257221,  0.00810255,  0.04392574,  0.09202785]],\n\n       [[ 0.2251054 ,  0.06344208, -0.05926704,  0.04094871],\n        [ 0.03012719,  0.14752444,  0.09756479,  0.03783219],\n        [-0.01962293,  0.04596839,  0.16733417,  0.02646091],\n        [ 0.01209812,  0.00708544,  0.04848414,  0.0750044 ]],\n\n       [[ 0.20561136,  0.07068361, -0.05526942,  0.04214355],\n        [ 0.02634501,  0.11640777,  0.08698693,  0.03201595],\n        [-0.02096767,  0.03766458,  0.14711928,  0.02239741],\n        [ 0.01240291,  0.01207825,  0.0407007 ,  0.06205465]]])"
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Phis"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[0.69406084, 0.19866282, 0.07542932, 0.12525187],\n        [0.19866282, 0.88773112, 0.06679637, 0.19003002],\n        [0.07542932, 0.06679637, 1.42715106, 0.06422977],\n        [0.12525187, 0.19003002, 0.06422977, 0.72200096]],\n\n       [[0.71726523, 0.21086386, 0.07762947, 0.13556694],\n        [0.21086386, 0.92522197, 0.09989765, 0.19786279],\n        [0.07762947, 0.09989765, 1.48927624, 0.07490291],\n        [0.13556694, 0.19786279, 0.07490291, 0.73020439]],\n\n       [[0.78261253, 0.23370429, 0.07995612, 0.15431813],\n        [0.23370429, 0.97427757, 0.12628964, 0.2116062 ],\n        [0.07995612, 0.12628964, 1.5816593 , 0.09332389],\n        [0.15431813, 0.2116062 , 0.09332389, 0.76725946]],\n\n       [[0.85317743, 0.2698687 , 0.06926143, 0.17331434],\n        [0.2698687 , 1.02196557, 0.13101231, 0.23299451],\n        [0.06926143, 0.13101231, 1.69132926, 0.10093577],\n        [0.17331434, 0.23299451, 0.10093577, 0.80573988]],\n\n       [[0.94045189, 0.30923633, 0.07197224, 0.1972989 ],\n        [0.30923633, 1.11821147, 0.19374696, 0.26337299],\n        [0.07197224, 0.19374696, 1.83387395, 0.13730869],\n        [0.1972989 , 0.26337299, 0.13730869, 0.84210924]],\n\n       [[0.98444108, 0.33271292, 0.06952747, 0.2118588 ],\n        [0.33271292, 1.15520034, 0.22149974, 0.27822534],\n        [0.06952747, 0.22149974, 1.90125637, 0.15162861],\n        [0.2118588 , 0.27822534, 0.15162861, 0.85656501]],\n\n       [[1.03418221, 0.35740251, 0.06540288, 0.22665428],\n        [0.35740251, 1.19070286, 0.24605777, 0.2937144 ],\n        [0.06540288, 0.24605777, 1.96600415, 0.16601709],\n        [0.22665428, 0.2937144 , 0.16601709, 0.87223308]],\n\n       [[1.07844559, 0.38103029, 0.06125979, 0.23996286],\n        [0.38103029, 1.22009087, 0.26634066, 0.30774115],\n        [0.06125979, 0.26634066, 2.02290186, 0.17944297],\n        [0.23996286, 0.30774115, 0.17944297, 0.88457293]],\n\n       [[1.11781577, 0.40192754, 0.059341  , 0.25223498],\n        [0.40192754, 1.25087273, 0.29181306, 0.32196036],\n        [0.059341  , 0.29181306, 2.0759561 , 0.19364168],\n        [0.25223498, 0.32196036, 0.19364168, 0.89508382]],\n\n       [[1.15067023, 0.41998803, 0.0566537 , 0.26267098],\n        [0.41998803, 1.273563  , 0.30987657, 0.3332263 ],\n        [0.0566537 , 0.30987657, 2.1170029 , 0.20392287],\n        [0.26267098, 0.3332263 , 0.20392287, 0.90308569]]])"
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.mse(10)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[0.69406084, 0.19866282, 0.07542932, 0.12525187],\n        [0.19866282, 0.88773112, 0.06679637, 0.19003002],\n        [0.07542932, 0.06679637, 1.42715106, 0.06422977],\n        [0.12525187, 0.19003002, 0.06422977, 0.72200096]],\n\n       [[0.73390433, 0.21019034, 0.08078239, 0.13415639],\n        [0.21019034, 0.92486466, 0.10191765, 0.1996849 ],\n        [0.08078239, 0.10191765, 1.48022615, 0.08068433],\n        [0.13415639, 0.1996849 , 0.08068433, 0.73037717]],\n\n       [[0.8137135 , 0.23007019, 0.08717984, 0.15714781],\n        [0.23007019, 0.9613197 , 0.119603  , 0.21042024],\n        [0.08717984, 0.119603  , 1.57101013, 0.10930129],\n        [0.15714781, 0.21042024, 0.10930129, 0.77287167]],\n\n       [[0.90254594, 0.26189279, 0.06732058, 0.17878345],\n        [0.26189279, 1.00363239, 0.13902531, 0.22821483],\n        [0.06732058, 0.13902531, 1.67790352, 0.11255464],\n        [0.17878345, 0.22821483, 0.11255464, 0.80922474]],\n\n       [[0.99691438, 0.29995471, 0.05751131, 0.1900392 ],\n        [0.29995471, 1.12097784, 0.2125487 , 0.25832551],\n        [0.05751131, 0.2125487 , 1.80684022, 0.1540173 ],\n        [0.1900392 , 0.25832551, 0.1540173 , 0.84472027]],\n\n       [[1.05884017, 0.3191863 , 0.05230346, 0.20161467],\n        [0.3191863 , 1.15975476, 0.24735837, 0.27256036],\n        [0.05230346, 0.24735837, 1.86623974, 0.17078755],\n        [0.20161467, 0.27256036, 0.17078755, 0.85804679]],\n\n       [[1.12268523, 0.34070397, 0.0444961 , 0.2121616 ],\n        [0.34070397, 1.1968314 , 0.27743765, 0.28799341],\n        [0.0444961 , 0.27743765, 1.92452431, 0.18883691],\n        [0.2121616 , 0.28799341, 0.18883691, 0.87334247]],\n\n       [[1.1817892 , 0.36024958, 0.03494831, 0.21919167],\n        [0.36024958, 1.23284467, 0.3065218 , 0.30153704],\n        [0.03494831, 0.3065218 , 1.97311774, 0.20351061],\n        [0.21919167, 0.30153704, 0.20351061, 0.88364266]],\n\n       [[1.2328996 , 0.37739745, 0.0280261 , 0.22399154],\n        [0.37739745, 1.27442603, 0.3396899 , 0.31571953],\n        [0.0280261 , 0.3396899 , 2.01680251, 0.21847345],\n        [0.22399154, 0.31571953, 0.21847345, 0.89227024]],\n\n       [[1.27886066, 0.39136703, 0.02203036, 0.22907553],\n        [0.39136703, 1.30337859, 0.36435224, 0.32671168],\n        [0.02203036, 0.36435224, 2.05020594, 0.22959503],\n        [0.22907553, 0.32671168, 0.22959503, 0.89865438]]])"
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "denominator.cumsum(axis=0)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "$\\mathrm{MSE}(h) = \\sum_{i=0}^{h-1} \\Phi \\Sigma_u \\Phi^T$"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0.69406084, 0.88773112, 1.42715106, 0.72200096],\n       [0.71726523, 0.92522197, 1.48927624, 0.73020439],\n       [0.78261253, 0.97427757, 1.5816593 , 0.76725946],\n       [0.85317743, 1.02196557, 1.69132926, 0.80573988],\n       [0.94045189, 1.11821147, 1.83387395, 0.84210924],\n       [0.98444108, 1.15520034, 1.90125637, 0.85656501],\n       [1.03418221, 1.19070286, 1.96600415, 0.87223308],\n       [1.07844559, 1.22009087, 2.02290186, 0.88457293],\n       [1.11781577, 1.25087273, 2.0759561 , 0.89508382],\n       [1.15067023, 1.273563  , 2.1170029 , 0.90308569]])"
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mse"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "outputs": [],
   "source": [
    "newde = np.zeros((10,4))\n",
    "for i in range(10):\n",
    "     newde[i] = np.diag(denominator[i])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[0.69406084, 0.88773112, 1.42715106, 0.72200096],\n       [0.73390433, 0.92486466, 1.48022615, 0.73037717],\n       [0.8137135 , 0.9613197 , 1.57101013, 0.77287167],\n       [0.90254594, 1.00363239, 1.67790352, 0.80922474],\n       [0.99691438, 1.12097784, 1.80684022, 0.84472027],\n       [1.05884017, 1.15975476, 1.86623974, 0.85804679],\n       [1.12268523, 1.1968314 , 1.92452431, 0.87334247],\n       [1.1817892 , 1.23284467, 1.97311774, 0.88364266],\n       [1.2328996 , 1.27442603, 2.01680251, 0.89227024],\n       [1.27886066, 1.30337859, 2.05020594, 0.89865438]])"
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newde.cumsum(0)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "results.orth_ma_rep()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[0.69406084, 0.04445819, 0.00398667, 0.02172854],\n        [0.05686377, 0.88773112, 0.00312634, 0.05001573],\n        [0.00819753, 0.00502602, 1.42715106, 0.00571393],\n        [0.02260325, 0.04067832, 0.0028907 , 0.72200096]],\n\n       [[0.72531671, 0.05752855, 0.00405574, 0.03085299],\n        [0.05748081, 0.91091438, 0.01847822, 0.0505951 ],\n        [0.00909142, 0.00949646, 1.47691608, 0.00787395],\n        [0.02377623, 0.04136701, 0.00687122, 0.72653139]],\n\n       [[0.80057737, 0.07184923, 0.00417043, 0.04047302],\n        [0.05922299, 0.9463773 , 0.01965633, 0.051259  ],\n        [0.0096836 , 0.01292839, 1.56577829, 0.00946249],\n        [0.02626311, 0.04591285, 0.01336346, 0.76391844]],\n\n       [[0.87784968, 0.09396406, 0.00594993, 0.04764221],\n        [0.0644747 , 0.98472746, 0.02421604, 0.05634608],\n        [0.00970776, 0.01297705, 1.6721782 , 0.00993667],\n        [0.02855078, 0.0504489 , 0.01347985, 0.79963248]],\n\n       [[0.96504998, 0.10700517, 0.00786715, 0.05501127],\n        [0.07780115, 1.08183498, 0.04749856, 0.07054381],\n        [0.01083662, 0.01996242, 1.79634859, 0.01380236],\n        [0.03033597, 0.05152778, 0.02383049, 0.82653923]],\n\n       [[1.01601017, 0.12413713, 0.00926922, 0.0637574 ],\n        [0.08125794, 1.10841577, 0.0615578 , 0.07389697],\n        [0.01110607, 0.02294636, 1.85358151, 0.0155555 ],\n        [0.03170686, 0.05325709, 0.02714888, 0.83682099]],\n\n       [[1.070173  , 0.13965937, 0.01096824, 0.07115581],\n        [0.08591954, 1.13615895, 0.07187451, 0.07837104],\n        [0.01124944, 0.02553635, 1.90995233, 0.01704994],\n        [0.03304568, 0.05516023, 0.03128386, 0.84842642]],\n\n       [[1.11795494, 0.15549975, 0.01349403, 0.07789   ],\n        [0.09098433, 1.1597257 , 0.08469827, 0.08365443],\n        [0.0113247 , 0.02735934, 1.95706374, 0.01839425],\n        [0.03395887, 0.05618167, 0.03472602, 0.85568245]],\n\n       [[1.1593488 , 0.16788649, 0.01607877, 0.0839249 ],\n        [0.09657926, 1.18529251, 0.10127205, 0.08957647],\n        [0.01151355, 0.03053733, 1.99813363, 0.02020227],\n        [0.03471154, 0.05695405, 0.03870934, 0.86069008]],\n\n       [[1.19525249, 0.18100423, 0.01827148, 0.08996964],\n        [0.10047253, 1.20164013, 0.11422502, 0.09363373],\n        [0.01158081, 0.03265457, 2.02972514, 0.02146143],\n        [0.03539943, 0.05781819, 0.04156274, 0.86433031]]])"
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numerators.cumsum(axis=0)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[0.69406084, 0.04445819, 0.00398667, 0.02172854],\n        [0.05686377, 0.88773112, 0.00312634, 0.05001573],\n        [0.00819753, 0.00502602, 1.42715106, 0.00571393],\n        [0.02260325, 0.04067832, 0.0028907 , 0.72200096]],\n\n       [[0.72531671, 0.05752855, 0.00405574, 0.03085299],\n        [0.05748081, 0.91091438, 0.01847822, 0.0505951 ],\n        [0.00909142, 0.00949646, 1.47691608, 0.00787395],\n        [0.02377623, 0.04136701, 0.00687122, 0.72653139]],\n\n       [[0.80057737, 0.07184923, 0.00417043, 0.04047302],\n        [0.05922299, 0.9463773 , 0.01965633, 0.051259  ],\n        [0.0096836 , 0.01292839, 1.56577829, 0.00946249],\n        [0.02626311, 0.04591285, 0.01336346, 0.76391844]],\n\n       [[0.87784968, 0.09396406, 0.00594993, 0.04764221],\n        [0.0644747 , 0.98472746, 0.02421604, 0.05634608],\n        [0.00970776, 0.01297705, 1.6721782 , 0.00993667],\n        [0.02855078, 0.0504489 , 0.01347985, 0.79963248]],\n\n       [[0.96504998, 0.10700517, 0.00786715, 0.05501127],\n        [0.07780115, 1.08183498, 0.04749856, 0.07054381],\n        [0.01083662, 0.01996242, 1.79634859, 0.01380236],\n        [0.03033597, 0.05152778, 0.02383049, 0.82653923]],\n\n       [[1.01601017, 0.12413713, 0.00926922, 0.0637574 ],\n        [0.08125794, 1.10841577, 0.0615578 , 0.07389697],\n        [0.01110607, 0.02294636, 1.85358151, 0.0155555 ],\n        [0.03170686, 0.05325709, 0.02714888, 0.83682099]],\n\n       [[1.070173  , 0.13965937, 0.01096824, 0.07115581],\n        [0.08591954, 1.13615895, 0.07187451, 0.07837104],\n        [0.01124944, 0.02553635, 1.90995233, 0.01704994],\n        [0.03304568, 0.05516023, 0.03128386, 0.84842642]],\n\n       [[1.11795494, 0.15549975, 0.01349403, 0.07789   ],\n        [0.09098433, 1.1597257 , 0.08469827, 0.08365443],\n        [0.0113247 , 0.02735934, 1.95706374, 0.01839425],\n        [0.03395887, 0.05618167, 0.03472602, 0.85568245]],\n\n       [[1.1593488 , 0.16788649, 0.01607877, 0.0839249 ],\n        [0.09657926, 1.18529251, 0.10127205, 0.08957647],\n        [0.01151355, 0.03053733, 1.99813363, 0.02020227],\n        [0.03471154, 0.05695405, 0.03870934, 0.86069008]],\n\n       [[1.19525249, 0.18100423, 0.01827148, 0.08996964],\n        [0.10047253, 1.20164013, 0.11422502, 0.09363373],\n        [0.01158081, 0.03265457, 2.02972514, 0.02146143],\n        [0.03539943, 0.05781819, 0.04156274, 0.86433031]]])"
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "irfs"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[ 0.00000000e+00, -2.77555756e-17,  2.77555756e-17,\n         5.55111512e-17],\n       [ 0.00000000e+00,  1.11022302e-16,  2.77555756e-17,\n        -8.32667268e-17],\n       [ 1.38777878e-17,  2.77555756e-17,  0.00000000e+00,\n         1.38777878e-17],\n       [ 5.55111512e-17, -8.32667268e-17,  2.08166817e-17,\n         1.11022302e-16]])"
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CovUU / deviations - sigma_u / sd_u"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[[ 0.83310314,  0.21085111,  0.06314009,  0.14740605],\n        [ 0.23846126,  0.94219484,  0.05591365,  0.22364197],\n        [ 0.09054019,  0.07089443,  1.19463428,  0.07559055],\n        [ 0.15034377,  0.20168866,  0.05376522,  0.84970639]],\n\n       [[ 0.1767933 ,  0.11432567, -0.00831089,  0.09552196],\n        [ 0.02484018,  0.15226048,  0.12390271,  0.02406999],\n        [ 0.02989801,  0.06686132,  0.22308075,  0.04647597],\n        [ 0.03424887,  0.02624291,  0.06309135,  0.06730848]],\n\n       [[ 0.27433678,  0.11966904,  0.01070915,  0.09808174],\n        [ 0.04173945,  0.18831602,  0.03432372,  0.02576635],\n        [ 0.02433476,  0.05858269,  0.29809765,  0.03985659],\n        [ 0.04986857,  0.06742289,  0.08057447,  0.19335731]],\n\n       [[ 0.27797896,  0.14871056, -0.0421841 ,  0.08467108],\n        [ 0.07246866,  0.19583196,  0.06752559,  0.07132374],\n        [ 0.00491565,  0.00697618,  0.32618999,  0.02177561],\n        [ 0.04782965,  0.0673502 ,  0.01078833,  0.18898158]],\n\n       [[ 0.29529698,  0.11419767, -0.04378614,  0.08584324],\n        [ 0.11544024,  0.3116208 ,  0.1525861 ,  0.11915423],\n        [ 0.03359849,  0.0835785 ,  0.35237819,  0.06217468],\n        [ 0.04225151,  0.03284634,  0.10173809,  0.16403279]],\n\n       [[ 0.22574364,  0.1308891 , -0.03744423,  0.09352073],\n        [ 0.05879447,  0.16303618,  0.11857167,  0.05790649],\n        [ 0.01641489,  0.05462541,  0.23923403,  0.04187051],\n        [ 0.03702555,  0.04158494,  0.05760544,  0.10139898]],\n\n       [[ 0.23272907,  0.12458829, -0.04121915,  0.08601401],\n        [ 0.06827595,  0.16656283,  0.1015712 ,  0.06688848],\n        [ 0.01197388,  0.050892  ,  0.2374254 ,  0.03865799],\n        [ 0.03658986,  0.04362495,  0.06430383,  0.10772854]],\n\n       [[ 0.21859081,  0.1258586 , -0.05025723,  0.08206216],\n        [ 0.07116731,  0.15351466,  0.11324207,  0.07268694],\n        [ 0.00867508,  0.04269653,  0.21705161,  0.03666476],\n        [ 0.03021895,  0.03196011,  0.05866995,  0.08518232]],\n\n       [[ 0.20345482,  0.11129571, -0.05084037,  0.07768456],\n        [ 0.07479929,  0.15989625,  0.12873917,  0.07695482],\n        [ 0.01374245,  0.05637359,  0.20265708,  0.04252081],\n        [ 0.02743479,  0.02779156,  0.06311349,  0.0707646 ]],\n\n       [[ 0.18948269,  0.11453269, -0.04682639,  0.07774798],\n        [ 0.06239606,  0.12785781,  0.11381112,  0.06369662],\n        [ 0.00820083,  0.04601356,  0.17774   ,  0.03548467],\n        [ 0.02622768,  0.0293964 ,  0.0534173 ,  0.06033429]],\n\n       [[ 0.18042366,  0.10836217, -0.04706061,  0.07302376],\n        [ 0.06239179,  0.12077659,  0.10908415,  0.06357924],\n        [ 0.00699124,  0.04424565,  0.16464642,  0.03395905],\n        [ 0.02364847,  0.02673619,  0.0531013 ,  0.05290264]]])"
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orth_irfs"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[1.16079718],\n       [1.15850561],\n       [1.02205437],\n       [1.11178523]])"
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fe.sum(1)[:, None]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[1.16079718],\n       [1.15850561],\n       [1.02205437],\n       [1.11178523]])"
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[ 0.18948269,  0.11453269, -0.04682639,  0.07774798],\n       [ 0.06239606,  0.12785781,  0.11381112,  0.06369662],\n       [ 0.00820083,  0.04601356,  0.17774   ,  0.03548467],\n       [ 0.02622768,  0.0293964 ,  0.0534173 ,  0.06033429]])"
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(Phis[h] @ CovUU) / deviations"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "outputs": [
    {
     "data": {
      "text/plain": "array([[ 0.18948269,  0.11453269, -0.04682639,  0.07774798],\n       [ 0.06239606,  0.12785781,  0.11381112,  0.06369662],\n       [ 0.00820083,  0.04601356,  0.17774   ,  0.03548467],\n       [ 0.02622768,  0.0293964 ,  0.0534173 ,  0.06033429]])"
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Phis[h] @ (CovUU / deviations)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}